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ff^ ■ Abstract. We calculate curvature perturbations in the scenario in which the 

curvaton field decays into another scalar field via parametric resonance. As a result 
of a nonlinear stage at the end of the resonance, standard perturbative calculation 
. techniques fail in this case. Instead, we use lattice field theory simulations and the 

I separate universe approximation to calculate the curvature perturbation as a nonlinear 

function of the curvaton field. For the parameters tested, the generated perturbations 
O ^, are highly non-Gaussian and not well approximated by the usual /nl parameterisation. 

^ ' Resonant decay plays an important role in the curvaton scenario and can have a 

%^ I substantial effect on the resulting perturbations. 
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PACS numbers: 98.80. Cq, 11.15.Kc 
1. Introduction 
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! Observations of tlie cosmic microwave background radiation are consistent with 
Gaussian perturbations [1], but there are tantalising hints of non-Gaussianity [2] at 



0\ ', a level that would be clearly observable with the Planck satellite and other future 
experiments |3] . This would rule out the simplest inflationary models with slowly rolling 
, , , scalar fields, which can only produce very-near ly- Gaussian perturbations |3H7]. 
^ ! Models which can generate large non-Gaussianities either during, or at the end 

', of, inflation include those with non-canonical Lagrangians P-IT6]. those where there is 
breakdown of slow-roll dynamics during inflation due to sharp features in the potential 
pTtlTS]. those with multiple fields with specific inflationary trajectories in the field 
space [TU1 - E7] and those where additional light scalars affect the dynamics of perturbative 
or non-perturbative inflaton decay [28ti39] . 

A well-known example of a multi-field model where large non-Gaussianities can 
be generated after the end of inflation is the curvaton scenario [lOHll]. In this model, 
the primordial perturbations arise from the curvaton field, a scalar field which is light 
relative to the Hubble rate during inflation but, unlike the inflaton, gives a subdominant 
contribution to the energy density. 
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After inflation the relative curvaton contribution to the total energy density 
increases. This affects the expansion of space, and the curvaton perturbations become 
imprinted on the metric fluctuations. The standard adiabatic hot big bang era is 
recovered when the curvaton eventually decays and thermalises with the existing 
radiation. The mechanism can be seen as a conversion of initial isocurvature 
perturbations into adiabatic curvature perturbations during the post-inflationary epoch. 
If the curvaton remains subdominant even at the time of its decay, its perturbations 
need to be relatively large to yield the observed amplitude of primordial perturbations. 
Consequently, in this limit the curvaton scenario typically generates significant non- 
Gaussianities |I5] - I52] . 

Almost all the analysis of the curvaton scenario is based on the assumption of a 
perturbative curvaton decay. As pointed out in refs. [53l|5l], it is, however, possible 
that the curvaton decays through a non-perturbative process analogous to inflationary 
preheating This is a natural outcome in models where the curvaton is coupled 

to other scalar fields which acquire effective masses proportional to the value of the 
curvaton field. 

During this short period of non-equilibrium physics, part of the curvaton condensate 
decays rapidly into quanta of other light scalar fields. In analogy to infiationary 
preheating, the decay is not complete and the remaining curvaton particles have to decay 
perturbatively. The properties of primordial perturbations generated in the curvaton 
scenario depend sensitively on the dynamics after the end of infiation until the time 
of the curvaton decay [581460] . As the preheating dynamics are highly nonlinear, it is 
natural to ask if this stage could generate large non-Gaussianities. 

Our aim in this work is to calculate the amplitude and non-Gaussianity of 
perturbations from curvaton preheating. To do this we use the separate universe 
approximation [ETHM] and classical lattice field theory simulations [551 ES]- This 
method was applied to inflationary preheating in refs. [37H39] . We solve coupled fleld 
and Friedmann equations to determine the expansion of the universe during the non- 
equilibrium dynamics and the fraction of the curvaton energy which is converted into 
radiation. From these, we calculate the curvature perturbation as a function of the 
local value of the curvaton fleld. We flnd that, at least for our choice of parameters, this 
dependence is highly nonlinear, implying very high levels of non-Gaussianity. 

The structure of the paper is as follows. In section [2] we review the curvaton 
model and derive useful general results. In section |3] we present an analytic calculation 
of production of curvature perturbations using linearised approximation of preheating 
dynamics, and demonstrate that the answers depend on quantities that are not 
calculable in linear theory. A full nonlinear computation is, therefore, needed. In section 
Hlwe show how this can be done, and compute the amplitude and non-Gaussianity of 
the perturbations for three sets of parameters. Finally, we present our conclusions in 
section 
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2. General theory 

2.1. Model 

We study the curvature perturbations generated in a simple curvaton model with 
canonical kinetic terms and the potential 

Vi^, a, x) = V{<P) + \m'a' + \g^a\' , (1) 

where is the inflaton, a is the curvaton and x is a scalar field which the curvaton will 
decay into. During infiation the infiaton potential V{(f)) dominates the energy density. 
We assume standard slow roll infiation, so the infiaton field rolls down its potential until 
it reaches a critical value 0* at which the slow roll conditions fail and infiation ends. 
As usual, the perturbations of the infiaton field generate curvature perturbations, but 
we assume that their amplitude is negligible. This requires ^ 10~^, where H is 

the Hubble rate and e is the slow roll parameter, both of which are determined by the 
infiaton potential 1^(0). 

Instead, we assume that the observed primordial perturbations are generated solely 
by the curvaton field cr, which should be light during infiation so that it develops nearly 
scale-invariant quantum fiuctuations similar to those of the infiaton field. The curvaton 
field remains nearly frozen at the value a* it had at the end of infiation, until the time 
tosc at which the Hubble parameter has decreased to if ~ m and the curvaton starts 
to oscillate around the minimum of its potential. These fiuctuations contribute to the 
energy density, and, therefore, they affect the expansion of the universe and generate 
curvature perturbations. 

Eventually, the curvaton a decays into lighter degrees of freedom. In the standard 
curvaton scenario this decay is assumed to take place perturbatively. The interaction 
term in ([1]) only allows pair annihilations, the rate of which falls quickly when the 
universe expands. Some of the curvaton particles would, therefore, survive until today, 
behaving as dark matter. While this might be an interesting scenario to study, we follow 
most of the literature and assume that the curvaton is also coupled to other fields. In 
particular, a Yukawa coupling to a light (possibly Standard Model) fermion field i\} of 
the form 

-^^Yukawa = ho^^ (2) 

would give the decay rate [ST] 

Without further knowledge of the curvaton couplings, the decay rate F appears as a 
free phenomenological parameter, which is bounded from below because the curvaton 
cannot decay arbitrarily late without spoiling the success of the hot big bang cosmology. 
The most conservative lower bound for the decay time is given by nucleosynthesis: 

r > t|bn/Mpi [S 
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In this work, we assume that the other scalar field x iii © is heavy during 
inflation, ga > H, so that its value remains close to zero and it develops no long- 
range perturbations. It was recently shown in ref. [53] that in this case some of the 
curvaton particles can decay via a parametric resonance into quanta of the x field. 
The resonance does not destroy all the curvaton particles, so a Yukawa coupling is still 
needed to allow those which remain to decay perturbatively. However, it reduces their 
density significantly. 



2.2. Curvature perturbations 



Using the 6N approach, the curvature perturbation ( on superhorizon scales is given by 
the expression [62] 



( = 6lna 



(4) 



where a|p is the scale factor at some fixed energy density p, and 5 denotes the difference 
from the mean value over the whole currently observable universe. The scale factor a is 
normalised to be constant at some earlier flat hypersurface, which we choose to be at 
the end of inflation. This measures the difference in the integrated expansion between 
Friedmann-Robertson-Walker (FRW) solutions with different initial conditions which 
are evolved until the same flnal energy density p. In our case, the variations of initial 
conditions arise from the super-horizon fluctuations of the curvaton fleld (5ct^, produced 
during inflation. This means that the curvature perturbation is a local function of the 
curvaton field fluctuations ( = ({6(J^). As we know the statistics of 6a^, this determines 
the statistics of the curvature perturbation ( completely. In this paper we compute this 
function. 

Observations show that the primordial perturbations were nearly Gaussian, which 
implies that C{So^^,) should be approximately linear. Therefore, it is natural to Taylor 
expand (jl]) as 

^ ^ ' (5) 



C = (lna)' 



Sa^ + ^(Ina)" 



where the prime denotes a partial derivative with respect to the curvaton value at the 
end of inflation a^, evaluated at constant flnal energy density p. This implies that, to 
leading order in 5cr*, the power spectrum of the curvature perturbations is 

2 



(In a)' 



(6) 



where Va{k) is the power spectrum of the curvaton fleld. Assuming the curvaton 
perturbations at the end of inflation 5a^, are Gaussian, the expansion ^ is of the 
form |i68j 

3 

5' 



c 



Cg + ^/nlc| + • • • 



(7) 
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where Q is a Gaussian field and /nl is a position independent constant. This defines 
the local nonlinearity parameter /nl and, according to (jS]), the 6N formalism gives 

5 (In a)" 



NL 



6 (In a)'2 



(8) 



Neglecting the inherent non-Gaussianity of the curvaton perturbations Sa^ yields an 
error proportional to slow-roll parameters. This is irrelevant in the limit |/nl| ^ 
which is the main focus of this work. The WMAP data and other recent observations 
give the constraint |/nl| < 100 [IH3l[70]. 

In this paper we focus on the case in which the curvaton field decays into x particles 
through a parametric resonance, at the end of which the fields undergo a period of 
potentially very complicated, nonlinear, non-equilibrium dynamics. Afterwards, the 
fields equilibrate, and we assume that eventually the universe behaves as a mixture of 
non-interacting matter and radiation. This assumption is checked using simulations in 
section H] and we find it to be accurate enough for our current purposes. Therefore, we 
parameterise the total energy density as a sum of matter and radiation components 



P = Pref 



if) 



+ (1 



(9) 



where Oref, Pref and rref = Pm,ref/Pref are the scale factor, energy density and matter 
fraction at some arbitrary reference time well after the resonance is over||| 

We assume that the infiaton field has either decayed into ultrarelativistic degrees of 
freedom or is itself ultrarelativistic, so that it contributes to the radiation component. 
The X fidd is also ultrarelativistic; the only degree of freedom that contributes to the 
matter component is the curvaton. 

If the resonance destroyed all of the curvaton particles, rref would be zero, but, in 
practice, some curvatons are left over. We assume rref ^ 1 which corresponds to the 
curvaton being subdominant at the end of the resonance. We can then rearrange (|9]) to 
give, at leading order in rref, the scale factor at energy density p 

1/4 



In a = In a^ef + 



\n — + r,ei 
P 




(10) 



where we have assumed that the curvaton is still subdominant at this time. That is 

1/4 



'^ref 



Pref 
P 



< 1. 



:iii 



The curvature perturbation is given by combining with (fTOj) . In general, Oref, Pref 
and rref all vary between one separate universe and another. In our case, the variation 
is entirely due to fiuctuations of a*, the value of the curvaton field during infiation. In 
order to calculate the curvature perturbation using we differentiate ([TU]) with respect 



I We define r = Pm/p following ref. 143 . This differs from the variable r used in ref. 
of 3/4 in the limit r <C 1. 



by a factor 
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to cr* keeping the final energy density p fixed. This gives 



(In a)' 



(In ttref )' + ^ 



(In Oref )' + - 



Prei I Pref 

ref , n I A 
Pref 4pV4pJ/4 

Pref 



1/4 



(In a) 




(12) 




(In a 



refj 



+ 



4 



Pref 
Pref 



/ \ 2 
Pref \ _|_ I Pref 

Pref, 
'ref 



4 \ Pref 




, . V ^refj-- . (13) 

'"ref Pref '"ref 

The remaining curvatons decay within a fixed time l/F where the decay rate F 
depends on the curvaton interactions. We assume that all matter in the universe 
is ultrarelativistic after the curvaton decay, so that the universe becomes radiation 
dominated and evolves adiabatically. Following the sudden decay approximation |15] . 
we assume that this decay takes place instantaneously at a fixed energy density pdecay, 
which is determined by the curvaton decay rate F. The final curvature perturbation is 



therefore given by setting p = pdecay in (fT2|) and f[T3|l . or equivalently r 

1/4 

_ / Hrci \ 

'"decay '"ref 



'"d 



ecay 5 



f Pref A 
\ Pdecay / 



where 



(14) 



The decay rate F is unknown, so we can treat rdecay as a free parameter. 



2. 3. Standard perturbative decay 

As a simple example of the calculation of perturbations, and to aid comparison with 
our results for resonant curvaton decay, we first summarise the standard perturbative 
calculation of perturbations generated in the curvaton model [l3|[i5] . 

This calculation ignores the resonance, so is assumed to be valid from the 
start of the curvaton oscillations, which we can therefore choose as our reference time, 
i'^ref = ^osc- This time is determined by the condition H ^ m, which implies that 



m'^Mpy independently of cr*. Furthermore, the scale factor a^ef 



IS 



Pref Post 

also independent of a* to good approximation, because the curvaton contribution to 
the energy density is negligible during inflation. Perturbations are therefore generated 
solely by derivatives of r^ef = rose, the matter fraction at the start of the oscillations, 
and (IT^ and (IT^ simplify to 



(In a)' 



'"decay '", 



'y osc 
'"nsr 



(15) 



Non-Gaussianity from resonant curvaton decay 7 

(In a)" = , (16) 

p 4 rose 

where we have also assumed rose ^ '"decay ^ 1- Using the results f l24|) and fl25|) below, 
the derivatives of the matter fraction rose = "^^CTosc/^Pose read 
r' 2 r" 2 

^ = _ , ^ = _ (17) 

^osc ^ose '^it: 

giving the simple result |i45j 
5 1 



/i 



NL 



3 '"deeay 

3. Linearised calculations 
3.1. Resonance 

Our general result in f|T2|) and f lT3l) depends on the quantities Oref, Pref and Vj-ef and 
their derivatives. Ideally, we would like to be able to calculate them analytically for 
general parameter values. We first attempt to carry out this calculation in linear theory 
of resonance [57], i.e. neglecting the backreaction of particles produced during the 
resonance. This approximation gives a clear physical picture of the resonance, and also 
an accurate quantitative description of many aspects of it. However, we find that it is 
not suitable for calculating the curvature perturbation, and, therefore, this calculation 
acts, ultimately, as a motivation for the nonlinear approach in section HI 

As we assume the x field is heavy during inflation and has a vanishing vacuum 
expectation value, and that the universe becomes radiation dominated after the end of 
inflation, the equation of motion for the curvaton is given by 
3 

o- + — cr + m^cr = , (19) 

which can easily be put into the canonical form of a Bessel equation. The general 
solution of ( !T9|) which remains bounded as t — )■ is given by 

a(t) = 2'/«r(5/4).. 1^ , (20) 

where Ji/4(x) is a Bessel function of the first kind and a{0) = cr* denotes the initial 
curvaton value set by the dynamics during infiation. In the model ([T]) that we consider 
here, the value of the curvaton field is constrained by 
H^, HlMl 
m 

where the upper bound comes from the subdominance of the curvaton and the lower 
bound is required to keep the x field massive during infiation and to enable broad 
parametric resonance 

The curvaton remains nearly frozen at the value cr* until the Hubble parameter 
decreases to H m and the field starts oscillate around the minimum of its potential. 
We assume the curvaton is still subdominant at this time and obeys the solution fl20|) . 
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After the onset of oscillations, mt > 1, fl20l) can be approximated by the asymptotic 
expression 

Following the notation of [57], we choose to normalise the scale factor to unity at 
^osc = 37r/ (8m) 



tosc J V Svr 

Up to corrections 0{H/m), the energy density of the oscillating curvaton is given by 
2 a 

where cXosc is the envelope of the oscillatory solution f l22|) at t, 



(23) 



_ ^OSC _ FCT.OSC /r)A\ 

Pf^ ~ r, „3 ~ ^3 ' ^ ^ 



^osc= y^l^ff a, ^0.76a, ■ (25) 

The time tosc appears as an unphysical reference point in our analysis but as it formally 
corresponds to a quarter of the first oscillation cycle it can also be thought of as the 
beginning of curvaton oscillations, hence the subscript. 

The coherent curvaton oscillations induce a time-varying effective mass for the x 
field which can lead to copious production of x particles as discussed in The 
process is analogous to the standard inflationary preheating scenario (SHIETIITT] , except 
that the universe is dominated by radiation instead of non-relativistic matter during the 
resonance. The resonant curvaton decay is efficient if 

g = ^»l, (26) 

which corresponds to broad resonance bands in momentum space. For curvaton values 
in the range 0211) this condition is always satisfied, as one can immediately see by taking 
into account the lightness of the curvaton during infiation, m -C if*. 

Using the solution fl22]) and applying the analytic methods developed to analyse 
particle production during the parametric resonance [57], the comoving number density 
of X particles at late times (t ^ tosc) can be estimated as 

nJt) ^ -—$■^=6^""'' , (27) 

where fc^, = ((^mcTosc)^^^ and /i ~ 0.14 is an effective growth index. Due to the stochastic 
nature of parametric resonance in expanding space, the actual growth index depends 
very sensitively on other parameters and varies rapidly between and 0.28 between one 
cycle of oscillation and another. 

As the number density of x particles grows exponentially, the contribution to the 
effective curvaton mass g'^{x^) eventually comes to dominate over the bare mass m?. 
The time tbr when the two contributions are equal — after which the backreaction of 
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the X particles can no longer be neglected — is found by equating the number density 
% ~ {x'^)g<y with ni^a/g. Using the result f l27l) . this yields an equation for tbr 

tb. - ^In ( IQV-/-^^-)^^^ , , (28) 



4m/i \ g^a, 



osc 



whose solution is given by the —1 branch of Lambert W-function 

tb. - -^w^-i (-10- W/^g^.e) • (29) 



Here gosc denotes the value of the resonance parameter fl26l) at the beginning of 
oscillations. The result is consistent with the assumption t ^ tosc for all yU provided 
that we place the mild constraint g^^^qUc < 10, which corresponds to tbr ^ lOtosc- In 
this limit (129 p can be approximated to reasonable accuracy by 

tbr ~ ^In (I0i°(7~ VCsc) + O (m-i In(mtbr)) ■ (30) 

The value of the scale factor at the time of backreaction is found by substituting 
m\i into ^ 

a.r - (-^) W'_{' (-10- VV/^g^4^) • (31) 

The derivatives of In a^r with respect to cr*, which are needed for computing the curvature 
perturbation, are 

(Inabr)' ^ - i(ln/.)' - (l + 0{a^^')) (32) 

(Inab.)" ^ - ^(In/.)" + ^^^^^ (l + 0{a^^')) , (33) 

where ' = d/da^ and the relation 9"/9cr" = {aosc/ d"^ / da^^^ following from (125|) has 
been used. The terms involving derivatives of the effective growth index /i would be 
very difficult to estimate analytically and we therefore leave them in an implicit form. 

Approximate energy conservation |m^(T^j, + 'T^^c^r — 1/2 "^^o'osc'^br^ gives an 
estimate 

« = (34) 

12 m^ag^ 3 ag^. 

for the resonance parameter (l26l) at the onset of backreaction. If Q'br ^ 1, the resonance 
does not terminate immediately at tbr- Instead, the effective curvaton mass m^^ = g'^{x^) 
starts to grow exponentially, which soon brings the resonance to end. The nonlinear 
stage of the resonance after tbr therefore makes only a small contribution to the total 
duration of the resonance [57] which can be reasonably well estimated by the duration 
of the hnear stage tbr given by (12^ . 
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3.2. Curvature perturbation 



It is not possible to describe the backreaction and the equihbration of the fields using 
linear theory. To carry out the linearised calculation, we therefore simply assume 
that the resonance ends instantaneously at time tbr, and that some fraction ^ of the 
curvaton energy is transferred into radiation. The matter and radiation energy densities 
immediately after backreaction are therefore given by 

Pa,osc 



Pm,hr 



Pr,br 



(1-0^ 

Pr,osc 



''br 



(35) 

(36) 



''br 



In order to calculate the curvature perturbation, we choose this backreaction time 



as our reference time in 

Pref = Phr " 
'"ref '"br " 



1]). Therefore we have Orcf = Obr, and 

Pm,hr + Pr .br 



Pr. osc , Per, osc 
— '■ \- ■ ^ 



Po 



(1 + rosc(«br - 1)) , (37) 



''br 



''br 



Pm, br 



(1 0Po-,oscObr 



'1 - OWflbr 



(38) 



Pm, br ~l~ Pr, br Pr, osc ~l~ Pa-,osc'-''br 

Pcr,osc/Pr,osc and we have assumed that roscObr ^ 1, which is 



where rose = P<7,osc/Posc 
equivalent to having r ^ 1 throughout the resonance. Substituting these into (IT^ and 
f|T3|l . we can work out the expression for the curvature perturbation. Using the leading 
terms in fl5^ . (1551) . we find 

2 e' 



(In a)' 
(In a)" 



^ Qbr 
^ Qbr 

1-e a 



2p e 



+ 



(39) 



2 4 
— + — 

0"* 



at 



27i^ 



+ V + 



p 



-1/2 



Cp 



2 



4^' 



0"* 



1-0 1-^ ■ 

By setting ^ = in the above expressions, the standard curvaton results (IT5]) for a 
perturbative decay are recovered. 

The terms scaling as in fl39l) and fHOl) represent the contribution due to radiation 
inhomogeneities created by the curvaton decay. If the curvaton particles left over after 
the resonance have a long perturbative lifetime, this component effectively redshifts 
away and the expressions (155]) and (110|) reduce to 



(In a)' 
(In a)" 



r 

r 
2^ 



^4' 



2(1-0 



2 c// 



l-e 2(1-0 



The nonlinearity parameter /nl can be computed using 
we find using (gT]) and , 



(41) 
(42) 

In the late time limit, 
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To make use of these results, one would have to know the first and second derivatives 
of fi and ^, which, unfortunately, are difficult to calculate. In principle, the exponential 
growth rate n is calculable in linearised theory, but obtaining its derivatives reliably is 
hard because it is an averaged quantity that describes a stochastic process. In contrast, ^ 
is fully determined by the nonlinear dynamics at the end of the resonance, and, therefore, 
it is not possible to calculate it using linearised equations. 

We conclude that the perturbations generated by resonant curvaton decay are not 
calculable using linear theory. Instead, we have to carry out a fully nonlinear calculation 
using lattice field theory methods. 



4. Simulations 

4.1. Simulation method 

In section E] we saw that the curvature perturbations produced by the resonance 
cannot be calculated using linear theory. Therefore, we adopt a completely different 
approach. We use nonlinear three-dimensional classical lattice field theory simulation^ 
to determine how the scale factor evolves as the fields fall out of equilibrium at the 
end of preheating and to compute the required quantities Oref, Pref and rref directly. 
This is a standard method of solving the field dynamics in such systems [SSlinSlE2HZZ] , 
and recently [371439] it was shown how it can be used, with the separate universe 
approximation, to compute curvature perturbations. 

The a and x fields are taken to be position-dependent, and the scale factor a is 
homogeneous over the whole lattice. The latter assumption is justified when the lattice 
is smaller than the Hubble volume, and allows us to describe the expansion of the 
universe using the Friedmann equation. We have a coupled system of field equations, 

a + 3H& - 4 W + mV + g^ax^ = , (44) 
a"' 

X + 3Hx-^V\ + g'a\ =0. (45) 

and the Friedmann equation 

r2 _ P_ 

where the energy density is calculated as the average energy density in the simulation 
box, 



The coupled system of equations (1441) . fl4S]) and (jlHD are solved on a comoving lattice 
with periodic boundary conditions. The infiaton component is included as idealised 
radiation where p^o is the density in the field at the beginning of the simulation. We 



§ We use a modification of the corrected version of the code used in refs. [37] and [3S] described in the 
erratum of 1551. 
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(a) a field (b) x field 



Figure 1. Evolution of the fields during one simulation run. To the left of the vertical 
dashed line the evolution is calculated with ((6T|) - (l65]) : to the right it is calculated with 
(|48)) - (|52)) . The cj) field is included as a homogeneous radiation component. 



solve the system in conformal time {dr = a ^dt) using a fourth order Runge-Kutta 
algorithm. The system of equations is then 



a" + 2-a' - — VV + (mV + g^ax^) 



X +2— X 
a 







. 



The second order Friedmann equation is 

a" = l^a^ 
6 ' 

where the averaged density and pressure are 
P 



1 

3a4 + Zs 



^ [icr'r + ixr) + ^ ((Va)^ + (Vx)^) + Via, x) 



d^x 



L2a2 



{^'f + (X ) 



/n2 



1 



(48) 
(49) 

(50) 

(51) 
(52) 



The initial infiaton energy density p^Q depends on the infiaton potential V{(f)). For 
simplicity, we assume it has a quartic form, 

1 



V(0) = ^A0^ 



(53) 



We choose the initial value of cf) to be (^q = Mpi, where we use the subscript to denote 
the initial values for the simulation. We still have the freedom to choose the value of 
the conformal time at the start of our simulations. We take 
^Mpi 



= V Alf • (^^' 

SO that the scale factor is simply proportional to conformal time, a(r) oc r even at 



T ~ To- 
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The initial values for the a and x fields consist of a homogeneous background 
component ctq or Xo^ which represents fiuctuations with wavelength longer than the 
lattice size, and Gaussian inhomogeneous fiuctuations which represent short-distance 
quantum fiuctuations. 

The X field is massive during infiation; its long-wavelength fiuctuations are heavily 
suppressed. In contrast, the curvaton field a is still frozen at the start of the simulation. 
As a result of superhorizon fiuctuations, its local value a* at the end of infiation varies 
between one Hubble patch and another; this sets the initial value o"o for the curvaton 
field in our simulation, (Tq = cr^. By running simulations with different initial curvaton 
values ctq, we can, therefore, determine how a^ei, Pref and Vj-ei depend on ctq and calculate 
the curvature perturbation using (jS]). 

The relevant range of ctq is determined by the power spectrum Vaik) of the curvaton 
field. Assuming that the curvaton mass m is very small, the power spectrum is given 
by 

where Hk and Nk are the Hubble rate and the number of e-foldings before the end 
of infiation evaluated at the time when the mode k = aH^ left the horizon. A more 
accurate calculation for non-zero m is given in appendix 1 of [38J. 

To completely determine the observable curvature perturbation, we need to run the 
simulation for all values of a* that were realised in our current Hubble volume. The 
curvaton field is a Gaussian random field, so the mean value over a Hubble volume 
today, dZ, has a Gaussian distribution. The variance of this distribution is 

rNtot / 1 \ 4 

^^'^ ^ Jn V ~ N'J "^^^ ^^^Pi^tot • (56) 

Note that this depends on A^tot? the total number of e-foldings of infiation, so is essentially 
a free variable. This distribution is centred around zero, assuming no homogeneous 
classical curvaton background. 

The width of the range of cr* that we need to cover is given by the variance of 
cr* within our current Hubble volume, between different Hubble volumes at the end of 
infiation. This is given by 

{S<tI) = j^' V. (k) (^1 - i-^ dN, ^ ^,XMl,Nl , (57) 

where Nq ^ 60 is the number of e-foldings after the largest currently observable scales 
left the horizon. This gives the width of the range of a^, that we need to consider. Note 
that the only free parameter this depends on is A. 

For the initial values of the inhomogeneous field modes, we follow the standard 
approach p^66p72p73p75] . The x field is given random initial conditions from a Gaussian 
distribution whose two-point functions are the same as those in the tree-level quantum 
vacuum state, 

(XkXq) = (2vr)^5(k-q)-^ , 
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(XkXq) = (27rmk-q)^, (58) 

where = ■\/^^~+~mJ = + q'^ctq- All other two-point correlators vanish. The 
inhomogeneous modes of a are populated similarly to those of x- 

In order to estimate the error bars in, for example, figure 3 we repeat each 
simulation with the same parameters multiple times using different realisations of the 
above random initial conditions. In earlier work on preheating [57H5^ . the whole 
evolution was calculated using full nonlinear equations. In the current case that would 
be computationally very expensive, because the model is not conformally invariant and 
the universe expands by a large factor ~ 100 during the evolution. To fit the relevant 
wavelengths inside the lattice throughout the whole simulation, the lattice size would 
have to be much larger than 100^, which is not possible given the number of simulations 
which need to be run. 

Instead, we take a shortcut, and make use of the fact that until the magnitude of the 
fluctuations in x become large enough to backreact on a (when (x^) ~ m'^/g'^) the field 
dynamics are linear to a very good approximation. Therefore, we can evolve the whole 
probability distributions of the initial conditions in fc-space using linearised versions of 
( H8|) -( l52|) . More precisely, we need to solve the linear equations for the inhomogeneous 
modes, 

a' 

a'i + 2—al + fcVk + a^mVk = , (59) 
a 

Xk + 2-xL + ^'Xk + a2^?V2xk = , (60) 
a 

in the background provided by the solutions of the nonlinear equations for the 
homogeneous modes, 

a" + 2%' + (mV + g'^ax^) = , (61) 

X" + 2-x' + a^g^a\ = , (62) 

a 

// Ip-BP 3 
" = 6^" ' 

P=^ + ^((0' + (xf)+n^,x), (64) 



The background solution is independent of the modes in fl59|l and fl60jl . and therefore it 
needs to be calculated only once for each (Tq. 

Crucially, the general solution of the linear mode ( 159|) and (160|) is obtained by 
solving the equations for two different initial values. This is because, as a consequence 
of linearity, the general solution can be represented as 

M(r) I Jjoi I ■ 
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where M(r) is a two-by-two matrix. Consider now the solutions with initial conditions 
of position 1 and velocity and vice versa: 



Xk(i,o)(^) 

Xk(l,0)(^) 




(67) 
(68) 



Xk(o,i)l^J 

Xk(0,l)(^) 

Linearity implies that the solution of arbitrary initial conditions Xk(0) and Xk(0) is given 
simply by 




-^<Hs;:::w)^^^'ns::::;w)- 

We use this approach to evolve the fields until some scale factor ai, which is well 
before the nonlinear terms become important. Only then, we Fourier transform the 
fields to coordinate space and start to evolve them using the full nonlinear equations 
(I48|) - (l52|) . This is shown in figure [H We checked that our results do not depend on the 
chosen value of ai. 

This method has two advantages which make the calculation presented here 
possible. Firstly, we evolve the whole distribution, so the linear evolution only needs 
to be calculated once, and only the relatively short nonlinear stage has to be repeated 
for each realisation of the initial conditions. Secondly, and more significantly, it reduces 
the factor by which the universe grows during the evolution of fl48p - fl52p from ~100 to 
~10. This allow us to use lattices with much lower resolution to achieve comparably 
accurate results, thus significantly reducing computation time. 



4-2. Simulation results 

Our model has three free parameters. A, m, g. It is not possible for us to probe this whole 
parameter space using current computational technology. Therefore, we simply pick one 
set of parameters, A = 10"^^ g = 0.01 and m = lO^^Mpi. We take xo = 10"^^Mpi. 

We also have to choose the initial value of the curvaton field cr. In order to cover 
all values of a* in our currently observable universe, we carried out simulations for a 
range of ctq 

oo - ^ScTo < (To < 0^ + ^Sao , (70) 
whose width is given by (1571) . 

ScTo ^ y^Wfi - lO-'Mpi . (71) 

The probability distribution of the central value Oq is given by f l56|) . but as a results of its 
dependence on the total number of e- foldings iVtot, which is unknown, we can essentially 
choose it freely. Nevertheless, the value of oq is restricted by several constraints. The 
resonance being brought to an end by backreaction demands that g 3> 1 (see ( 12^ ). 
The characteristic frequency of the x fi^ld oscillations is gao and the time the system 
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Figure 2. The matter fraction r during one simulation. The resonance is rapid and 
destroys ^ 95% of the matter. The two solid lines show the two definitions of r given 
in ([75)1 and ([71| . pref = 5 x 10"^'^ is marked by the vertical dashed line. The blue 
dashed line represents the approximation Q taken about this prof 

takes to reach backreaction is proportional to m~^. Therefore, the simulation time is 
proportional to gao/m = 2g^/^. In our simulations we used three values: 

= 0.0005Mpi, O.OOlMpi and 0.002Mpi . (72) 

These corresponds to qbr ~ 50, 200 and 750 respectively. With these parameters one 
32^ point lattice simulation takes ~3 hours. We repeat each simulation several (~50) 
times with different random realisations of initial conditions to obtain statistical errors. 
Although, this computational method is very resource intensive, the independence 
of each simulation (or Hubble patch) means that the problem scales perfectly on 
a multiprocessor machine. A time-step of (2 x 10~^/m)Mp^ and lattice spacing of 
{5m/96)Mp^^ were used. 

Each simulation gives us time-streams of the scale factor a, the density p, given by 
f lSTj) . and the pressure P, given by f l52|) . We measure the matter fraction r = pm/ P in 
two different ways: from the ratio of the pressure and energy density, 

P 

r = 1 - 3- , (73) 
P 

and by fitting the dependence of the scale factor on the energy density, 

log a = —- log p + constant, (74) 

in a narrow range about the point of interest. In continuous time these two expressions 
would agree, and therefore any discrepancy between them would be a sign of time 
discretisation errors. 
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Figure 3. Results from simulations for different ranges of ctq centred around 
ctq — O.OOOSAfpi, CTo = O.OOlMpi and ctq = 0.002Mpi, each covering the values present 
in one Hubble volume today. The top and bottom rows show In Orcf and Tj-ef respectively 
measured at p = 5 x 10~^'^Mpj and averaged over 10 — 50 runs. The solid lines show 
quadratic fits of the form (I75p . The centre panels also show some numerical checks of 
the simulation results. The green (square) points represent simulations with time-steps 
double the length of the black points. The blue (diamond) points are represent a four 
times longer time-step. It can be seen that the results are within the errors of the 
primary runs. 



The measured matter fraction r is shown in figure [2] as a function of the energy 
density. The two curves correspond to the two definitions fl73|) and fl7i|) . which clearly 
agree very well. Initially, it grows because the matter component decreases more slowly 
than the radiation component. At ^ 10^°Mpi^, it drops rapidly as the resonance 
destroys most of the curvatons. Afterwards r again grows as the universe expands. 

To calculate the curvature perturbation using ( IT^ and (IT^ . we need to choose 
a reference point at which we measure the scale factor Oref, energy density pref and 
matter fraction r^^f. We chose this be at fixed energy density pref = 5 x 10"^^Mpj. We 
interpolated our measurements using a simple least-squares fit in {log(p), log(a)} and 
{log(p), log(r)} about p = prcf to find the scale factor a^ef and the matter fraction rref. 
The results are shown in figure [3] for three choices of oq. 
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We found that the choice of pref has no effect on the results as long as it is well after 
the end of the resonance. This is demonstrated by figure [21 in which the (blue) dashed 
curve shows the matter fraction r calculated from the assumption ([9]) of non-interacting 
matter and radiation components with the measured values of Oref and rj-ef. This agrees 
very well with the measured r{p) at late times. 

We calculate the first and second derivatives of InOref and r^^f with respect to ctq. 
Assuming that the expansion works well over the whole range of cxo in the current 
Hubble volume, we do this by fitting quadratic polynomials 



Inaref(o-o) = Inaref(cro) + \na[^f{ao - (^o) + ^^^a"A^o 
rrc{{(^o) = rrciio^) + r',^{{(^o - o^) + ^rfcf (o"o - o^)^ 



^)^ 



(75) 



to our data. These fits are shown by the curves in figure [3] and the corresponding fit 
parameters are given in table [TJ 

We concentrate on two measurable quantities: the spectrum of curvature 
perturbations Vt^, which is given by (|6]), and the nonlinearity parameter /nl, which 
is given by ([8]). We chose to fix pref to be a constant, so f|T2|) and f|T3l) simplify to 



(In a)' 
(In a)" 



(InOref)' + 



'"decay '"ref '"ref 

4 r^ef 



Hti ^ V _L '"fiecay ''ref '"ref 
(^maref j + - 



4 Tref 

For the amplitude of the curvature perturbations to agree with observations, ~ 10 



(76) 
(77) 



-10 



and (jTHD imply 



ecay 



r^ef + 4— 



ref 




Using this, the nonlinearity parameter (|8]) can be recast in the form 



/i 



NL 




— (In a 




(78) 



(79) 



The spectrum of the curvaton fluctuations produced during inflation is given by (155|) . 
which, for our parameters is, P^. ^ 5 x 10~^'^Mp[, so that \fVcJVl ~ 45. 

The values of rdecay and /nl calculated from the fit parameters are shown in table [2l 
In all cases rdecay has an acceptable value, r^ef <^ ?"dccay ^ 1, which means that the 
curvaton can produce perturbations with the observed amplitude. On the other hand, 
/nl is much larger than the observations allow [IHSIEO], which rules out these parameter 
values. 

Note that these values of rdecay and /nl rely on the quadratic fit fl75l) . which clearly 
does not describe the data very well, as can be seen in figure [31 This illustrates that the 
curvature perturbations produced in this model are not well approximated by truncating 
the expansion ([7|) at second order and parameterising the non-Gaussian effects only by 
local /nl- Therefore one should not put too much emphasis on the precise numbers, but 



Non-Gaussianity from resonant curvaton decay 19 



U 


0.0005 


0.001 


0.002 


(In ttref)' 


0.024 ±0.002 


-0.005 ±0.001 


0.025 ±0.016 


(Inflrcf)" 


57000 ± 13000 


-44000 ± 9000 


-420000 ± 120000 


'^ref 


(1.116 ±0.004) X 10"^ 


(3.72 ±0.002) X 10-6 


(10.4 ±0.011) X 10-^ 




0.109 ±0.008 


-0.009 ± 0.004 


-0.014 ±0.024 


r" 
' ref 


(2.1 ±0.6) X 10^ 


(1.7 ±0.3) X 10^ 


(-8.3 ±1.8) X 10^ 



Table 1. Fit parameters defined in ([75]) in Planck units. 





0.0005 


0.001 


0.002 


'"decay 


0.0018 ±0.0001 


0.08 ±0.03 


0.13 ±0.21 


/nl 


(3.6 ±1.0) X 10^ 


(-3.7 ± 1.7) X 10^ 


(-1.0 ± 1.8) X 10^ 


port 
'"decay 


0.046 


0.091 


0.182 


rpert 

/nl 


36 


18 


9 



Table 2. The values of rdccay and /nl calculated using (l78t and (l79t with the full 
simulation data (upper) and in the standard perturbative curvaton theory (lower). 



the conclusion that the predicted perturbations are highly non-Gaussian and ruled out 
by current observations. 

However, our calculations are non-perturbative and we are not limited to the usual 
Taylor expansion in ([5]). Instead, using ffTOj) we obtain the whole function C('^*) as a 
simple linear combination 

C((T*) = lna((j*) - In a(a-*) = (51naref(o"*) ± C(5ri.ef (cr*), (80) 

where the constant C depends on the perturbative curvaton decay rate T and is given 
by 



In figure H] we show this full result for the three choices of oq, calculated with 
the values of rdecay shown in Table [21 These results give a complete description of 
the statistics of curvature perturbations in this model, and they could be used to 
produce maps of the CMB or to compare directly with observations. The correct way to 
calculate /nl and other nonlinearity parameters describing higher order statistics from 
our results would, therefore, be to repeat the precise procedure the observers use in 
their measurements. The data are not well approximated by ([7j) truncated at second 
order, so different ways of doing this would probably yield very different values for /nl- 
Also, /nl and other nonhnearity parameters would also probably appear to be scale and 
position dependent, because on smaller scales one would cover a smaller range of cr*. 

Figure [3] shows that, with our parameters, Slnaj-^i is orders of magnitude smaller 
than the amplitude of the observed perturbations, and, therefore, the dominant 
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Figure 4. The curvature perturbation as a function of the curvaton field value for 
the three cases (from top to bottom) = O.OOOSMpi, O.OOlMpi, 0.002A/pi, calculated 
using ([50)1 with the value of rdccay taken from table [21 In each case, the dependence is 
highly nonlinear, implying strong non-Gaussianity. 



contribution to ([80]) must come from the second term, i.e. from the variation of the 
matter fraction Tref- Figure [2] shows that the matter fraction drops by a factor of around 
20 during the nonhnear stage, but in itself, this drop has no effect on the properties of 
the final perturbations, because it could be fully compensated for by reducing pdecay 
Instead, what is important is that the amount of curvatons left after the nonlinear stage 
depends sensitively on the value of a*. This leads to the variation of rref within the 
range of a^, present in one Hubble volume, which is seen in figure [3] and which makes 
the dominant contribution to the curvature perturbations. 

We believe that the sensitive dependence on a* is caused by the nonlinear field 
dynamics at the end of the resonance, but we do not have a detailed explanation 
of this. Nevertheless, it is in qualitative agreement with the results found in [571176] 
for inflationary preheating, and with findings that the resonance dynamics in massless 
preheating are effectively chaotic and that small variations in initial conditions can have 
a great impact on curvature perturbations [32] ■ In the current case the dependence 
seems smooth, suggesting that the dynamics are not strictly chaotic. 
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5. Conclusions 

In this work we have presented a method of calculating the curvature perturbations from 
resonant curvaton decay. The curvaton field is light relative to the Hubble rate during 
inflation and, as a result, its fluctuations are correlated on superhorizon scales; the 
mean value of the curvaton field varies between one Hubble volume and another. The 
local value of the curvaton field affects both the amount of expansion and the fraction 
of the curvaton particles which are destroyed during resonant decay into another field. 
This turns fluctuations of the curvaton field into curvature perturbations, which are 
generally non-Gaussian. As the field evolution at the end of the resonance is nonlinear, 
these perturbations cannot be calculated using standard perturbative techniques. 

We computed this effect numerically using classical field theory simulations 
combined with the 6N approach. This approach was used earlier in refs. [37H39] . but 
the current case is technically more difficult because the model is not scale invariant. 
To deal with this, we divided the evolution into three parts: First, the inhomogeneous 
modes were evolved numerically using linearised equations. Shortly before nonlinearities 
become important, we switched to a full nonlinear simulation to describe the non- 
equilibrium dynamics at the end of the resonance. Eventually, when the system has 
equilibrated, we extrapolated to late times by assuming that the universe consists of 
non-interacting matter and radiation. This approach allowed us to study evolution 
during which the universe expands by several orders of magnitude, which would not be 
possible using a straightforward field theory simulation. 

Our method is fully nonperturbative; we are not forced to assume that the non- 
Gaussianity is of the simple form ([7]) parameterised by /nl, and indeed we find that, at 
least at our parameter values, it is not. This may well be a fairly generic consequence of 
nonlinear field evolution. Instead, our method produces the full numerical dependence 
of the curvature perturbation ( on the curvaton field a*, which is a Gaussian random 
variable with a known power spectrum. More work is needed to properly compare such 
data with observations. 

Our findings demonstrate that the curvaton resonance can leave a significant 
imprint on primordial perturbations. At the parameter values we studied the 
perturbations were far too non-Gaussian to be compatible with observations, but it 
is very likely that there are other parameter values with acceptable levels of non- 
Gaussianity. Therefore it would be interesting to explore a wider range of parameters 
in future work. In this work we focused on parameters for which the curvaton is 
subdominant throughout the whole evolution and the decay product field x is massive 
during inflation, but the other possibilities are equally important and worth investigating 
in future work. We also expect a very similar resonance to take place in a model in 
which the curvaton is charged under a gauge group and decays resonantly into gauge 
bosons. 
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